# Marginal effects --------------------------------------------------------

interaction_plot_continuous_ggplot <- function(model, effect, moderator, interaction, varcov = "default", minimum = "min", maximum = "max", incr = "default", num_points = 50, conf = 0.95, mean = FALSE, median = FALSE, alph = 80, rugplot = TRUE, histogram = TRUE, title = "Marginal effects plot", xlabel = "Value of moderator", ylabel = "Estimated marginal coefficient") {
  
  # Extract Variance Covariance matrix
  if (varcov == "default") {
    covMat = vcov(model)
  } else {
    covMat = varcov
  }
  
  # Extract the data frame of the model
  mod_frame = model.frame(model)
  
  # Get coefficients of variables
  beta_1 = model$coefficients[[effect]]
  beta_3 = model$coefficients[[interaction]]
  
  # Set range of the moderator variable
  # Minimum
  if (minimum == "min") {
    min_val = min(mod_frame[[moderator]])
  } else {
    min_val = minimum
  }
  # Maximum
  if (maximum == "max") {
    max_val = max(mod_frame[[moderator]])
  } else {
    max_val = maximum
  }
  
  # Determine intervals between values of the moderator
  if (incr == "default") {
    increment = (max_val - min_val) / (num_points - 1)
  } else {
    increment = incr
  }
  
  # Create data frame of moderator values at which marginal effect is evaluated
  data <- data.frame(moderator = seq(from = min_val, to = max_val, by = increment))
  
  # Compute marginal effects
  data$delta_1 <- beta_1 + beta_3 * data$moderator
  
  # Compute variances
  data$var_1 <- covMat[effect, effect] + (data$moderator^2) * covMat[interaction, interaction] + 2 * data$moderator * covMat[effect, interaction]
  
  # Standard errors
  data$se_1 <- sqrt(data$var_1)
  
  # Upper and lower confidence bounds
  z_score = qnorm(1 - ((1 - conf) / 2))
  data$upper_bound <- data$delta_1 + z_score * data$se_1
  data$lower_bound <- data$delta_1 - z_score * data$se_1
  
  # Create ggplot
  p <- ggplot(data, aes(x = moderator, y = delta_1)) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    geom_line() +
    geom_ribbon(aes(ymin = lower_bound, ymax = upper_bound), alpha = 0.2) +
    labs(title = title, x = xlabel, y = ylabel)
  
  return (p)
}